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Abstract 



The numerical implementation of finite element discretization method 
.^ ■ for the stream function formulation of a linearized Navier-Stokes equa- 

^^ i tions is considered. Algorithm 1 is applied using Argyris element. Three 

global orderings of nodes are selected and registered in order to conclude 
the best banded structure of matrix and a fluid flow calculation is con- 
sidered to test a problem which has a known solution. Visualization 
of global node orderings, matrix sparsity patterns and stream function 
contours are displayed showing the main features of the flow. 
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1 Introduction 



The numerical treatment of nonlinear problems that arise in areas such as fluid 
mechanics often requires solving large systems of nonlinear equations. Many 
methods have been proposed that attempt to solve these systems efficiently; 
one such class of methods is the finite element method, the most widely used 
5^ I technique for engineering design and analysis. 

The Navier-Stokes equations may be solved using either the primitive vari- 
able or stream function formulation. Here, we use the stream function formu- 
lation. The attractions of the stream function formulation are that the incom- 
pressibility constraint is automatically satisfied, the pressure is not present in 
the weak form, and there is only one scalar unknown to be determined. The 
standard weak formulation of the stream function first appeared in 1979 in [9]; 
a general analysis of convergence for this formulation has been done in [4, 5]. 

The goal of this paper is to demonstrate that the method can be imple- 
mented to approximate solutions for incompressible viscous fiow problems. 



2 Governing Equations 



Consider the Navier-Stokes equations describing the flow of an incompressible 
fluid 



-Re^^ A M + (m- V)m+ Vp = f in n, 


(2.1) 


y -u = in fi, 


(2.2) 


■u = on dil, 


(2.3) 



curl/ in Q 


(2.4) 


on dn 


(2.5) 


on on, 


(2.6) 



where u = {ui,U2) and p denotes the unknown velocity and pressure field, 
respectively, in a bounded, simply connected polygonal domain fl C R"^] /is 
a given body force; and Re is the Renolds number. 

The introduction of a stream function ilj{x,y) defined by 

oy ox 

means that the continuity equation (2.2) is satisfied identically. The pressure 
may then be eliminated from (2.1) to give 

Re"^ a"^ iIj - iljy A ip^ + ip^ A ipy -- 

i, -. 
dip 
dh 

where h represents the outward unit normal to Vt. Equation (2.4) is a nonlinear 
fourth-order partial differential equation which turns into the known linear 
biharmonic fourth-order partial differential equation by omitting the second 
and third terms in the left-hand side of the equation. In order to write (2.4) - 
(2.6) in a variational form, we define the Sobolev spaces 

H^{Vl) = {v:veL'^{n),DveL'^{n)}, (2.7) 

H^{n) = {v:veH\n),v = on dVt}, (2.8) 

H\Q) = {v:veL\Q),Dve L\Q), D^v G L^}, (2.9) 

H^{n) = {v.ve H\n) ■.v = — = 0, on dQ}, (2.10) 

where L'^{Q) is the space of square integrable functions on Q and D represents 
differentiation with respect to x or y. For each ip G H^{Q), define curl (f = 

"^ J . The following linearized weak form of equations (2.4) - (2.6) is 

~'-Px J 

considered: 

Find ip G El{^) such that for all ^ G //o(fi), 
a(^,V9)+6(^*;^,<^)=%), (2.11) 



where 

a{^, if) = Re~^ / A ■?/' ■ A cpdn, 
Jn 

Jn 

i{(p) = (/, curl (p) = / / ■ curl (pdil, 

Jn 

where ip* is a fixed given function (a primitive approximation for ip). 

3 Finite Element Discretization 

For the standard finite element discretization of (2.11), we choose conforming 
finite element subspace X C Hq{Q). We then seek ip^ E X such that for all 

a(^^/)+6(^*^^^^^) = %'*). (3.1) 

The existence and uniqueness for the solution of the discrete problem (3.1) has 
been proved in [1]. Once the finite element spaces are prescribed, problem (3.1) 
reduces to solving a system of algebraic equations. Various iterative methods 
can be used to solve problem (3.1). In this paper the following algorithm 
has been applied to solve problem (3.1) for a fixed mesh of size h. In each 
iteration in Algorithm 1, we need to solve a linear system. The resulting linear 
system is nonsymmetric whose symmetric part is positive definite. Moreover, 
the resulting matrix is a sparse matrix. We choose the conjugate gradient 
stabilized (BIGG STAB) method, which requires two matrix- vector products 
and four inner products in each iteration. 



Algorithm 1 (the finite element algorithm) 

Given Max-iteration & Tolerance 
Given tp^ as a starting guess 
For i = 1 : Max- iteration 

Solve the linear system on the mesh X C H^ for ipi: 
a{tlJi,(p) + b{^lJi_i,tlJi,(p) =i{(p) \/ (f e X 

If \\ipi — ipi-i\\ < Tolerance & Residual < Tolerance 

Stop 
End 



4 Finite Element Space 

The inclusion X C Hq{Q) requires the use of finite element functions that 
are continuously differentiable over Q. We choose Argyris triangle as a finite 
element space for the stream function formulation. We will impose boundary 
conditions by setting all the degrees of freedom at the boundary nodes to 
be zero and the normal derivative equal to zero at all vertices and nodes on 
the boundary. In Argyris triangle, the functions are quintic polynomials within 



each triangle and the 21 degrees of freedom are chosen to be the function value, 
the first and second derivatives at the vertices, and the normal derivative at 
the midsides. 

5 Global Node Orderings 

For linear systems derived from a partial differential equation, each unknown 
corresponds to a node in the discretization mesh. Different orderings of un- 
knowns correspond to permutations of the coefficient matrix. Since the con- 
vergence speed of iterative method may depend on the ordering used, the next 
three options have been considered. 

1. Start numbering the nodes at the vertices with six successive numbers 
at each node. Then the numbering of the midpoint nodes is performed 
starting with the horizontal sides, followed by the vertical and oblique 
side. This is illustrated in Figure (5.1). 

2. Start numbering the function values of the nodes at vertices followed by 
the derivatives, and then go back to the midpoint nodes using the same 
ordering as above. This is illustrated in Figure (5.2). 

3. Start numbering the nodes at the vertices with six successive numbers by 
skipping a node alternatively, and then go back to the midpoint nodes of 
horizontal, vertical and oblique sides using the same orderings as above. 
This is illustrated in Figure (5.3). 
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Figure 5.2 
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Figure 5.3 



It can be concluded that the first option is the most appropriate since it 
has the best matrix property, being a banded structure, as shown in Figures 
5.4, 5.5, 5.6 which visuahze, respectively, the location of the nonzero elements 
for the three selected global orderings. Therefore this ordering has been used 
in the next section. 
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Figure 5.5 
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Furthermore, the following table illustrates the cpu-time, the number of 
arithmetic operations, and the number of bicgstab iterations which have been 
computed by performing the suggested global numberings in a specified prob- 
lem: 





cpu-time 


n.c.o. 


bicgstab iter. 


1. 


28.45 


227800974 


282.5 


2. 


29.77 


231648608 


289.0 


3. 


29.06 


239406016 


297.5 



6 Numerical Example 



In this section we describe some numerical results obtained by implementing 
the finite element algorithm, for which we have an exact solution. The region 
Q is the unit square {0<x<l, 0<|/<1} and for the finite element 
discretization, we use the Argyris elements. The biharmonic equation along 
with equation (2.11) was solved in the same way. They were solved until the 
norm of the difference in successive iterates and the norm of the residual were 
within a fixed tolerance. 

We consider as a test example the two-dimensional Navier-Stokes equations 
(2.1) - (2.3) on the unit square Q = (0, 1)^, where we define the right-hand side 
by / := — Re~^ A u + (u- V)u+ Vp with the following prescribed exact solution: 



u 






with il){x,y) = X {x — I) y {y — I) 



P 
7 



x^ + y^ -0.5. 



For this test problem, all requirements of the theory concerning the geometry 
of the domain and the smoothness of the data are satisfied. Moreover, the 
stream function ip{x, y) satisfies the boundary conditions of the stream function 
equation of the Navier-Stokes equations. 

In all numerical calculations in this example, we have used the Argyris ele- 
ments with Re = 1 and tol = 10~^. We pick three values of h: 1/3, 1/5, and 
1/9. The cpu-time, the number of PCG iterations, and the error {ipc — ipe\ for 
different values of h are tabulated in Table 6.1 based on four quadrature points, 
and Table 6.2 based on six quadrature points for the biharmonic equation. The 
cpu-time, the number of BICGSTAB iterations, errors, and the number of the 
mathematical operations are tabulated in Table 6.3 based on six quadrature 
points for the linearized equation (2.11), where ip* is computed by solving the 
biharmonic equation as initial guess, and where 

n.q.p. = number of quadrature points, 

n.c.o. = number of operations, 

pcg-itr. = number of PCG iterations, 

bicgstab itr. = number of BICGSTAB iterations, 

ipc = The computed solution of the function values, 

ipe = The exact solution of the function values. 



h 


n.q.p 


cpu time 


n.c.o. 


Error 


pcg-itr 


1/3 


4 


1.3200 


17601719 


0.1644 X 10-3 


72 


1/5 


4 


6.7500 


246112068 


0.1899 X 10-3 


211 


1/9 


4 


64.4300 


3.7401 X 10^ 


0.1376 X 10-3 


437 


Table 6.1 


h 


n.q.p 


cpu time 


n.c.o. 


Error 


pcg-itr 


1/3 


6 


1.8200 


16476551 


0.1984 X 10-3 


72 


1/5 


6 


7.4700 


249166334 


0.1935 X 10-3 


209 


1/9 


6 


64.7100 


3.6482 X 10^ 


0.1379 X 10-3 


454 



Table 6.2 



h 


cpu time 


n.c.o. 


^c-^eO 


\^c - ^e 1 


\^c - V'e 2 


bicgstab itr. 


1/3 


13.24 


17528466 


0.2589 X 10-3 


0.1294 X 10-2 


0.1692 X 10"^ 


100.5 


1/5 


35.32 


227800974 


0.2148 X 10-3 


0.1062 X 10-2 


0.1048 X 10-^ 


282.5 


1/9 


179.88 


3.50366282 x 10^ 


0.1423 X 10-3 


0.6986 X 10-3 


0.6016 X 10-2 


567.5 



Table 6.3 
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exact solution 



h = l/9 



Fig. 6.1: Streamlines for h = 1/3, 1/5, 1/9 with Re = 1 using finite element 
method on the linear problem with four quadrature points as shown in Table 

6.1 
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Fig. 6.2: Streamlines for h = 1/3, 1/5, 1/9 with Re = 1 using finite element 
method on the linear problem with six quadrature points as shown in Table 

6.2 
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Fig. 6.3: Streamlines for h = 1/3, 1/5, 1/9 with Re = 1 using finite element 
method on the linear problem with six quadrature points as shown in Table 

6.3 
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